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ABSTRACT 

A two dimensional oscillating flow analysis has been 
conducted simulating the gas flow inside Stirling engine heat 
exchangers. Both laminar and turbulent oscillating pipe flow 
has been investigat ed numerically for i?e max =1920 (Va = 80), 
10800 (Va 272), 19300 (Va = 272), 60800 (Va = 126). The 

results are compared with experimental results of previous 
investigators. Also predictions of the flow regime on present 
oscillating flow conditions have been checked by comparing 
velocity amplitudes and phase difference with those from lam- 
inar theory and quasi-steady profile. A high Reynolds number 
k — c turbulence model was used for turbulent oscillating pipe 
flow. Finally, performance evaluation of the k — e model was 
made to explore the applicability of quasi-steady turbulent 
models to unsteady oscillating flow analysis. 

NOMENCLATURE 

c fi — coefficient used in eddy viscosity equation 
c (i — turbulence model constants for e equation 
(a=U) 

D = diameter of the pipe, m 

k = turbulent kinetic energy, m 2 /s 2 , [= 1/2 (u n +r' 2 + te' 2 )] 
/ = turbulence length scale, m 

lm = mixing length, m, (=«y) 

L = length of the pipe, m 

P = pressure, N/m 2 

r = radial coordinate, m 

R = radius of the pipe, m 

Rtmaz = Reynolds number, (= |u m | ■ D/v) 

t = time, s 

u — axial velocity, m/s 

\u\ = amplitude of axial velocity, m/s 

v = velocity vector, m/s 

Va = Valensi number (dimensionless frequency), (= wR 2 /v) 
t = axial coordinate, m 

6 = Stokes-layer thickness, m/rad 1 ^ 2 , [= (2i//u>)ij 

€ = dissipation rate of turbulent kinetic energy, m 2 /s 3 

A = friction factor 


T} = dimensionless radial coordinate, (r/R) 
p = viscosity, Ns/m 2 
vj = kinematic viscosity, m 2 /s 
p = density, kg/m 3 

Gk — turbulent Prandtl number for k equation 
<r ( = turbulent Prandtl number for € equation 
r *= shear stress, N/m 2 
w angular frequency, rad/s 

Superscripts 

= time averaged value 
' = fluctuating velocity component 

Subscripts 
c — critical value 
cl — values at a pipe centerline 
i,j = denotes spatial coordinates [=(r,0,x)] 
m — cross-sectional mean value 
u’ = value at pipe wall 

INTRODUCTION 

The problems of periodic turbulent internal flow has 
been studied by many research workers experiment ally as well 
as computationally (I-£). The periodic flow can be divided 
into two classes: 1) unsteady flow with non- zero mean velocity 
and 2) unsteady flow with zero mean velocity. Previous re- 
searchers have called the first class of periodic flow, pulsating 
flow, and the second class, oscillating flow. 

The effect of imposed periodic pulsations on the time- 
averaged properties has been previously investigated. 
Ramaprian and Tu (1) found that for sufficiently high fre- 
quencies, the time-averaged flow variables, e g., velocity, wall 
shear stresses and power loss due to friction, were affected by 
the imposed unsteadiness. They also concluded that, from 
their computations, a quasi-steady turbulence modej cannot 
adequately describe unsteady flow conditions, at leaat for 
high frequencies. Ohmi et al. (3) performed extensive pul- 
sating flow experiments to determine transition criteria via 
non-dimensional parameters on the basis of velocity measure- 
ments. The early investigation of Sarpkaya (5), by detecting 
the growth rate of artificial disturbances, and Hershey ana 
Im (£), by friction factor measurements, show that critical 
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Reynolds number Re maI ,c is much different from the steady- 
flow value of about 2300. 

Hino et al. (7, §) investigated oscillating turbulent flow 
in pipes and rectangular ducts experimentally for various 
Reynolds numbers and frequencies. Many useful turbulent 
flow data have been generated through their work. They 
found that the turbulent energy production becomes exceed- 
ingly high in the decelerating phase, but the turbulence is 
reduced to a very low level at the end of the decelerating 
phase and in the accelerating stage of reversal flow. 

So far, no extensive numerical calculations have been 
given for oscillating flows covering wide ranges of Reynolds 
number and frequencies. Kohler (9) has performed numerical 
simulations of turbulent oscillating flows and compared mean 
velocity profiles and fluctuations with experimental results 
which were obtained from the oscillating now test facility at 
the University of Minnesota [see Seume and Simon (10)]. 

In the present paper, turbulent oscillating pipe flow has 
been analyzed by solving time- averaged continuity and mo- 
mentum equations using the standard k — t turbulence model 
of Launder and Sparding (11). A computer code, CAST, de- 
veloped by Peric and Scheuerer (12), has been modified to 
handle unsteady inlet conditions. Since most constants for 
this k — t model originated from steady turbulent flow mea- 
surements together with simplified assumptions, the validity 
of this quasi-steady turbulence model for unsteady flow con- 
ditions has been questioned. Performance evaluation of this 
k ~ e model for the present oscillating flow conditions will be 
discussed; One objective is to check the applicability of this 
quasi-steady turbulent model to oscillating flow analysis. 

ANALYSIS 

Governing Equations and Numerical Scheme 

The time-averaged equations of continuity and motion 
for an unsteady Newtonian fluid with constant, fluid proper- 
ties can be written as follows in vector form. 

Equation of continuity: 

(V ■ v) = 0 (1) 

Equation of motion: 

P^ = -V p -[VF] (2) 

In handling turbulent flow condition^. the above dependent 
variables represent time-averaged values and turbulent shear 
stress term in Equation (2) is closed by the “high Reynolds 
number version'’ of k-e turbulence model which requires so- 
lutions of the following transport equations: 

= ilit/ok) [V 2 *] +P-pt (3) 


P ^ [V 2 <] + - c (7 p<) (4) 


The present k-t model has the eddy viscosity given by 
= (5) 

and the production rate of turbulent kinetic energy is given 
by 


p £, < 6 > 

The empirical constants are the standard values, = 
0.09, c u = 1.44, c <2 = 1.92, <7* = 10, o ( = 1.3. The wall 
function method is used to specify wall boundary conditions. 

A conservative finite volume method with collocated 
variable arrangement by Peric et al. (1£) is used for dis- 
cretization. The convection fluxes in the model transport 
equations are discretized with the first-order upwind differ- 
encing scheme. Central differences are used to approximate 
the diffusion fluxes. 

Boundary and Initial Conditions 

Particular attention was paid to a flow geometry used 
by previous workers to perform the experimental oscillating 
flows; see Ohmi et al. (J). The flow geometry is a circular 
pipe, as shown in Figure (1), with diameter of 50.4 mm and 
length of 6000 mm; This gives LjD ^ 60 at the middle section 
of the pipe. At this middle section the flow can be considered 
fully developed. The inflow is taken to be uniform over the 
cross section and time dependent according to the relation: 

^m,in = j^ml Sinu?f (7) 


r 



Fig, 1 Geometry of the pipe and inlet condition 
for oscillating flow 


FLOW COMPUTATIONS AND RESULTS 

Four cases have been studied for 2000 < R( ma r < 60000 
with different oscillating frequencies [see Table (l)j. 62 x 
22 grid nodes were used in the axial and radial directions 
respectively. A grid independence test has been carried out 
by doubling the grid nodes. For each case, 180 time steps 
were used in one cycle and the results were collected after 4 
consecutive cycles. Convergence criteria has been set to 0.1% 
of the residual norms for every dependent variable and about 
2000 seconds of CPU time (Cray XMP) per cycle has been 
obtained. 

Ohmi et al. (14) derived the following critical Reynolds 
number, for fully aeveloped flow, up to which the theoretical 
correlation for laminar oscillating flow is valid. 

= 882 (8) 

VVa 

According to equation (8), Cases 1 and 2 are in laminar 
regime and Cases 3 and 4 fall into the turbulent regime [see 
Table (I)]. In the following sections, computational results for 
fully developed oscillating flow are analyzed for each test case 
by comparing with both experimental data (3, 14) and steady 
state analyses. 
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Table I An outline of test cases 


Testcase 


Va 

Re_/(V»r 

Eqn. (8) 

Re—.yo 7 *) 1 '* 

Eqn. (17) Hino el 

11.10) 

Case 1 

1920 

80 

215 

BBZ 

620 

77 

0 

Case 2 

10800 

272 

655 



676 



Case 3 

19300 

272 

1171 



676 



Case 4 

60800 

126 

5409 

V 

64 0 

\ 

/ 


Calculation of Oscillating Flow in Laminar Regime 

Figures (2a j and (2b) show the instantaneous velocity 
versus u )t at different radial locations for Cases 1 and 2 re- 
spectively. The plots in those figures are made for the axial 
location at the middle of the pipe. The plots show the in- 
crease in the velocity amplitude as the distance from the wall 
increases; also shown on the plot are differences in phase an- 
gle between the fluid motion near the wall and that near the 
pipe center line. Profiles of the velocity amplitude normalized 
by the center line value versus normalized radial distance for 
Case 1 and 2 axe shown in Figures (3a) and (3b). Also shown 
on the plots Eire the phase differences Detween the velocity at 
any raaial position and that of the center line. In those Fig- 
ures, Ohmi’s experimental data (3) for both amplitude ratio 
and phase differences are shown and an excellent agreement 
was obtained except, near the pipe wall large discrepancies in 
phase differences are noticed. It was found that, from laminar 



theory, phase difference increases rapidly as approaching the 
pipe wail and present numerical results follows the trend of 
laminar theory [see Figures (3a) and (3b)] and the same de- 
viation was found in earlier work between experimental and 
laminar theory [see (2)]. Consequently, the above discrep- 
ancies near the wall for phase differences can be attributed 
to difficulties in near-wall measurement. In Figure (3b), the 
dashed line represents steady turbulent velocity profile (1/7 
power law). This shows that Case 2, although on the bor- 
derline between laminar/turbulent regime, is very close to 
laminar flow condition and still away from turbulent motion 
[see Figure (6)]. 



V 

VJ 

I 

IT 

"si 


Fig. 3a Profile of the velocity amplitude and 
phase difference for Case 1 

(Velocity amplitude: 

solid line; present numerical data, 

AAA; experimental data, 

Phase difference: 

chain dot line; present, numerical data, 

-{- + +; experimental data.) 

Calculation of Oscillating Flow in Turbulent Regime 

Calculation were also made for oscillating flow in the 
turbulent regime. Inlet turbulent kinetic energy k tn was ob- 
tained assuming isotropic turbulence: 

fc.n = ^(77 X u m ,,n ) 2 ( 9 ) 


Fig. 2 Laminar axial velocity distribution ver- 
sus crank angle during a cycle 

(a) Case 1 

(b) Case 2 

(Data were obtained at r/R = 0, 0.25,0.55, 0.75, 
0.88, 0.97 from top to bottom respectively.) 


where TI is turbulent intensity at the inlet which is assumed 
to be 1% for Cases 3 and 4. Also the turbulence length iscale 
at the inlet was approximated by the following equation 
[see Peric and Scheuerer (12)]: 


lc 3/2 

hn = (10) 
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Fig. 3b Profile of the velocity amplitude and 
phase difference for Case 2 

(Velocity amplitude: 

solid line; present numerical data, 
dashed line; 1/7 power law profile, 
AAA; experimental data, 

Phase difference: 

chain dot line; present numerical data. 
+ + -f; experimental data.) 



0 6 12 0 6 12 

u ;t/(*r/6) u;*/(7r/6) 

Fig. 4 Turbulent axial velocity distribution ver- 
sus crank angle during a cycle 

(a) Case 3 

(b) Case 4 

(Data were obtained at r/R = 0, 0.25,0.55, 0.75, 
0.88, 0.97 from top to bottom respectively.) 


Equation (10) holds under local equilibrium conditions and 
for a logarithmic velocity law [see Rodi (15)] and in the equa- 
tion, the turbulent dissipation rate e in was estimated from 
Equation (5): 



(ID 


The following assumption for turbulent viscosity [Kohler (9)] 
vras used: 




^ R Cm a T C/i M 


( 12 ) 


Then Equation (11), after substituting Equation (12), be- 
comes 



With the above relationships, turbulence length scale was es- 
timated about 20% of the pipe diameter. 



Fig. 5a Profile of the velocity amplitude and 
phase difference for Case 3 

(Velocity amplitude: 

solid line; present numerical data, 
dashed line; laminar profile, 
dotted line; 1/7 power law profile, 
AAA; experimental data, 

Phase difference: 

chain dot line; present numerical data, 
•f 4- -f ; experimental data.) 
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t ig. 5b Profile of the velocity amplitude and 
phase difference for Case 4 

(Velocity amplitude: 

solid line; present numerical data, 
dashed line; 1/7 power law profile, 
AAA; experimented data, 

Phase difference: 

chain dot line; present numerical data, 
-f -f +; experimented data.) 


Computational results are shown in Figures (4) and (5) 
in the same order as in the laminar cases. Figures (4a) 
and (4b) show the time averaged velocity versus u it at differ- 
ent radial location for Cases 3 and 4 respectively. The plots 
show, in contrast with the laminar cases, much less variation 
in velocity phase with radial position. Specially, in Figure 4b, 
it is noticed that no phase differences are shown across the 
pipe cross section. Also profiles of the normalized velocity 
amplitude versus normalized radial distance for Cases 3 and 4 
are shown in Figures (5a) and (5b). Also shown on the plots 
are the phase difference Detween the velocity amplitudes. It 
should be noted that in Case 4 the normalized velocity was 
obtained using the amplitude of the “cross-sectional mean” 
velocity (rather than the center line one * used in Cases 1, 2, 
and 3). Also, the phase differences in Case 4 was obtained 
from the difference between the velocity at any radial posi- 
tion and the “cross-sectional mean” velocity (rather them the 
center line one - used in Case 1, 2, and 3). 

Experimental results for Case 3 indicate that there are 
turbulent bursts in decelerating phase near the wail (J). As 
seen in Figure (5a), numerical results for both velocity am- 
plitude ratio and velocity phase differences are different ap- 
preciably from the experimental data with maximum error of 
14% and 40%- respectively . Also, on Figure (5a), two ref- 
erence plots are shown for the normalized velocity, one for 
laminar flow (dashed line) and the other for 1/7 power law 


profile (dotted line). It can be seen that the velocity ampli- 
tude ratio obtained numerically (solid line) nearly matches 
with 1/7 power law profile. Experimental data show interme- 
diate value between the laminar and turbulent calculations in 
Figure (5a). Consequently, it indicates that the flow regime 
for Case 3 is in between laminar and turbulent (i.e., transient 
regime) and present high Reynolds number k - e turbulence 
model does not adequately predict the oscillating fluid motion 
in this regime. In Case 4, experimental results show' turbu- 
lent bursts on the velocity distributions during most of the 
cycle and agree well with steady turbulent correlation (14). 

In Figure (5b), numerical results also shows good agreement 
with experimental data (< 7%). Specially, in Figures (4b) 
and (5b), both computation and data agree that the velocity 
phase remains nearly constant. Consequently, it is concluded 
that the flow regime for Case 4 reaches the fully turbulent 
regime or turbulent quasi-steady state as mentioned by ear- 
lier researchers. 

CRITICAL REYNOLDS NUMBER 

Hino et al. (§) summarized the critical Reynolds number 
from several experimental and theoretical results. From his 
work, two critical Reynolds numbers were found for the oscil- 
lating pipe flow. Experimental critical Reynolds numbers are 
found for Res based on Stokes-layer thickness 6 = (2 
and the cross-sectional mean velocity amplitude lumj which 
can be expressed in terms of Re max and Va as in Equation (8): 

- V2Re 6 , c (14) 

v/Va v 

Using visualization and power measurement technique, 
Sergeev (16) has determined critical Reynolds number to be 
Res, c = 500 and also using hot wire measurement. Hino 
et al. (£) obtained Re 6fC = 550; these result in values for 
Re maj|C />/Ua equal 710 and 770 respectively. As mentioned 
earlier, Ohmi et al. (14) predicted a value of 882 [Equa- 
tion (8)1. They also derived another expression for critical 
Reynolds number using, instead of laminar oscillating flow 
theory, the following steady turbulent correlation together 
with the turbulence generation region assumption which was 
mentioned previously: 

T u . = \p\u m \ 2 /8 1 (15 ) 

A = 0.3164/Jfc^ J 

and obtained the following expression for critical Reynolds 
number: 

Re maIl c = (21lVVa) 8/7 (16) 


or 

Remote _ (2 n)g/7 y^r 1/7 (17) 

>/Va 

in terms of the expression of Equation (14). The above criti- 
cal Reynolds numbers are tabulated for each case in Table J. 

Table I shows the flow conditions for each case examined in 
this paper as w^ell as the estimated values of jyVa by 

different investigators. 

Also, plotted on Figure(6)are y/Va versus Re jnat show- 
ing the four cases examined in this study, as well as Equa- 
tions (8) and (17). From Figure (6), it is shown that Case 1 
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Fig. 6 Oscillating flow regime and critical 
Reynolds number. 


is laminar and Case 4 is turbulent according to available cor- 
relations [Eqn. (S), (IT)]. This finding is in agreement with 
the present numerical calculations. Case 2 is in the laminar 
flow regime according to Equation (8) and lies at the criti- 
cal line [Eqn. (IT)]: the present analysis shows that Case 2 is 
laminar. Case 3 is in the turbulent flow regime according to 
Equations (8) and (IT), however the current analysis shows 
that it lies in between laminar and turbulent. More extensive 
computation will be conducted to complete the map shown 
in Figure (6) from computational point of view as well as 
experimental data available today. 

SUMMARY AND CONCLUSIONS 

Numerical calculations have been performed for oscil- 
lating flows in laminar, transition and turbulent regimes and 
comparisons have been made with experimental data. For 
turbulent computations, Launder- Spalding k — e turbulence 
model has been employed with standard values of the model 
constants. The calculation method for the inlet values of 
turbulent kinetic energy k and turbulent length scale l has 
been discussed. Finally, several criteria for defining critical 
Reynolds number for oscillating pipe flow have been reviewed 
an3 compared according to the present numerical results. The 
conclusions for the results of the present investigation are 
given below; 

[1] Oscillating flow in laminar flow regime can be simulated 
numerically with relatively high accuracy. 

[2] Standard version of k - e turbulence model cannot ade- 
quately model oscillating flow in transition flow regime. 

[3] For fully turbulent regime or quasi- steady turbulent 
regime, the standard version of k - e turbulence model 
can predict the oscillating flow within an allowable er- 
ror. 


[4] Critical Reynolds number defined by Ohmi et al. (14) 
from his experimental results seems to be reliable one 
from numerical calculations, however, since definition of 
transition region is not clear and also present numerical 
test data are scarce, extensive computations are needed 
to establish more reliable criteria for Stirling engine ap- 
plications. 

[5] For more rigorous numerical calculation for the oscil- 
lating flow in transition regime which is the case for 
most Stirling Cycles, further improvement in turbulence 
modeling is also necessary. 
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